import math
import time

import numpy as np
import gurobipy as gp
from gurobipy import GRB

import verify_piecewise_utils
from NN import cdinn2
from prob import *

# case_name = "zhrR9"


def block_space(dim_min, dim_max):
    n = len(dim_min)

    sub_blocks = []

    # 递归地对当前维度进行划分
    def split_block(current_min, current_max, index):
        if index == n:
            sub_blocks.append([current_min, current_max])
            return
        mid = (current_min[index] + current_max[index]) / 2

        # 左子块
        left_max = current_max[:]
        left_max[index] = mid
        split_block(current_min, left_max, index + 1)

        # 右子块
        right_min = current_min[:]
        right_min[index] = mid
        split_block(right_min, current_max, index + 1)

    # 开始分割
    split_block(dim_min, dim_max, 0)
    return sub_blocks


def barrier_f_without_der(m, x0, wb_list, wb_length, useVars):
    x_layer = x0
    for i in range(wb_length):
        W = wb_list[i * 2]
        b = wb_list[i * 2 + 1]
        if useVars:
            y_layer = m.addMVar((W.shape[0],), vtype=GRB.CONTINUOUS, lb=-np.inf, ub=np.inf)
            E = np.identity(y_layer.shape[0])
            expr = W @ x_layer + b - E @ y_layer
            m.addConstr(expr == 0)
        else:
            y_layer = W @ x_layer + b

        if i != wb_length - 1:
            if useVars:
                z_layer = m.addMVar((y_layer.shape[0],), vtype=GRB.CONTINUOUS, lb=-np.inf, ub=np.inf)
                # z = activate(y)
                for j in range(y_layer.shape[0]):
                    m.addConstr(z_layer[j] == gp.max_(y_layer[j], 0.0))
            else:
                z_layer = y_layer
                for j in range(y_layer.shape[0]):
                    z_layer[j] = max(y_layer[j], 0.0)
            x_layer = z_layer
        else:
            x_layer = y_layer

        if useVars:
            m.update()
    return x_layer


def barrier_Bf_nonlinear(wb_list, wb_length, min_bound, max_bound):
    m = gp.Model()
    m.setParam('Outputflag', 0)
    m.setParam('NonConvex', 2)
    A = m.addVars(DIM, name="A")
    b_l = m.addVar(name="b_l")
    b = m.addVar(name="b")
    x0 = m.addMVar((DIM,), vtype=GRB.CONTINUOUS, name='x0', lb=-np.inf, ub=np.inf)
    for i in range(DIM):
        x0[i].lb = min_bound[i]
        x0[i].ub = max_bound[i]
    x1 = m.addMVar((DIM,), vtype=GRB.CONTINUOUS, name='x1', lb=-np.inf, ub=np.inf)
    for i in range(DIM):
        x1[i].lb = domain_min[i]
        x1[i].ub = domain_max[i]

    xx1 = m.addVar()
    xx3 = m.addVar()
    xx5 = m.addVar()
    xx7 = m.addVar()
    xx9 = m.addVar()
    xx11 = m.addVar()
    xx13 = m.addVar()
    xx15 = m.addVar()
    xx17 = m.addVar()
    m.addConstr(xx1 == x0[1] * x0[1])
    m.addConstr(xx3 == x0[3] * x0[3])
    m.addConstr(xx5 == x0[5] * x0[5])
    m.addConstr(xx7 == x0[7] * x0[7])
    m.addConstr(xx9 == x0[9] * x0[9])
    m.addConstr(xx11 == x0[11] * x0[11])
    m.addConstr(xx13 == x0[13] * x0[13])
    m.addConstr(xx15 == x0[15] * x0[15])
    m.addConstr(xx17 == x0[17] * x0[17])
    m.addConstr(x1[0] == x0[0] + time_step * ((x0[1] + 2 * x0[2] + 2 * x0[3] + 2 * x0[4] + x0[5]) / 100.0 + 1))
    m.addConstr(x1[1] == x0[1] + time_step * (x0[2]))
    m.addConstr(x1[2] == x0[2] + time_step * (-10 * (x0[1] - (xx1 * x0[1]) / 6.0) - x0[1]))
    m.addConstr(x1[3] == x0[3] + time_step * (x0[4]))
    m.addConstr(x1[4] == x0[4] + time_step * (-10 * (x0[3] - (xx3 * x0[3]) / 6.0) - x0[1]))
    m.addConstr(x1[5] == x0[5] + time_step * (x0[6]))
    m.addConstr(x1[6] == x0[6] + time_step * (-10 * (x0[5] - (xx5 * x0[5]) / 6.0) - x0[1]))
    m.addConstr(x1[7] == x0[7] + time_step * (x0[8]))
    m.addConstr(x1[8] == x0[8] + time_step * (-10 * (x0[7] - (xx7 * x0[7]) / 6.0) - x0[1]))
    m.addConstr(x1[9] == x0[9] + time_step * (x0[10]))
    m.addConstr(x1[10] == x0[10] + time_step * (-10 * (x0[9] - (xx9 * x0[9]) / 6.0) - x0[1]))
    m.addConstr(x1[11] == x0[11] + time_step * (x0[12]))
    m.addConstr(x1[12] == x0[12] + time_step * (-10 * (x0[11] - (xx11 * x0[11]) / 6.0) - x0[1]))
    m.addConstr(x1[13] == x0[13] + time_step * (x0[14]))
    m.addConstr(x1[14] == x0[14] + time_step * (-10 * (x0[13] - (xx13 * x0[13]) / 6.0) - x0[1]))
    m.addConstr(x1[15] == x0[15] + time_step * (x0[16]))
    m.addConstr(x1[16] == x0[16] + time_step * (-10 * (x0[15] - (xx15 * x0[15]) / 6.0) - x0[1]))
    m.addConstr(x1[17] == x0[17] + time_step * (x0[18]))
    m.addConstr(x1[18] == x0[18] + time_step * (-10 * (x0[17] - (xx17 * x0[17]) / 6.0) - x0[1]))

    Bf = barrier_f_without_der(m, x1, wb_list, wb_length, True)
    m.addConstr(A[0]*x1[0] + A[1]*x1[1] + A[2]*x1[2] + A[3]*x1[3] + A[4]*x1[4] + A[5]*x1[5] + A[6]*x1[6] + A[7]*x1[7] + A[8]*x1[8] + A[9]*x1[9] + A[10]*x1[10] + A[11]*x1[11] + A[12]*x1[12]
                + A[13]*x1[13] + A[14]*x1[14] + A[15]*x1[15] + A[16]*x1[16] + A[17]*x1[17] + A[18]*x1[18] + b_l <= Bf)
    m.addConstr(A[0]*x1[0] + A[1]*x1[1] + A[2]*x1[2] + A[3]*x1[3] + A[4]*x1[4] + A[5]*x1[5] + A[6]*x1[6] + A[7]*x1[7] + A[8]*x1[8] + A[9]*x1[9] + A[10]*x1[10] + A[11]*x1[11] + A[12]*x1[12]
                + A[13]*x1[13] + A[14]*x1[14] + A[15]*x1[15] + A[16]*x1[16] + A[17]*x1[17] + A[18]*x1[18] + b_l + b >= Bf)
    m.setObjective(b, GRB.MINIMIZE)
    m.optimize()
    if m.status == 2:
        return True, A[0].x, A[1].x, A[2].x, A[3].x, A[4].x, A[5].x, A[6].x, A[7].x, A[8].x, A[9].x, A[10].x, A[11].x, A[12].x, A[13].x, A[14].x, A[15].x, A[16].x, A[17].x, A[18].x, b_l.x + b.x
    else:
        return False, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None
def barrier_f_linear(wb_list, wb_length, min_bound, max_bound):
    m = gp.Model()
    m.setParam('Outputflag', 0)
    m.setParam('NonConvex', 2)
    A = m.addVars(DIM, name="A")
    b_l = m.addVar(name="b_l")
    b = m.addVar(name="b")
    x0 = m.addMVar((DIM,), vtype=GRB.CONTINUOUS, name='x0', lb=-np.inf, ub=np.inf)
    for i in range(DIM):
        x0[i].lb = min_bound[i]
        x0[i].ub = max_bound[i]
    B = barrier_f_without_der(m, x0, wb_list, wb_length, True)
    m.addConstr(A[0]*x0[0] + A[1]*x0[1] + A[2]*x0[2] + A[3]*x0[3] + A[4]*x0[4] + A[5]*x0[5] + A[6]*x0[6] + A[7]*x0[7] + A[8]*x0[8] + A[9]*x0[9] + A[10]*x0[10]
                + A[11]*x0[11] + A[12]*x0[12] + A[13]*x0[13] + A[14]*x0[14] + A[15]*x0[15] + A[16]*x0[16] + A[17]*x0[17] + A[18]*x0[18] + b_l <= B)
    m.addConstr(A[0]*x0[0] + A[1]*x0[1] + A[2]*x0[2] + A[3]*x0[3] + A[4]*x0[4] + A[5]*x0[5] + A[6]*x0[6] + A[7]*x0[7] + A[8]*x0[8] + A[9]*x0[9] + A[10]*x0[10]
                + A[11]*x0[11] + A[12]*x0[12] + A[13]*x0[13] + A[14]*x0[14] + A[15]*x0[15] + A[16]*x0[16] + A[17]*x0[17] + A[18]*x0[18] + b_l + b >= B)
    m.setObjective(b, GRB.MINIMIZE)
    m.optimize()
    if m.status == 2:
        return True, A[0].x, A[1].x, A[2].x, A[3].x, A[4].x, A[5].x, A[6].x, A[7].x, A[8].x, A[9].x, A[10].x, A[11].x, A[12].x, A[13].x, A[14].x, A[15].x, A[16].x, A[17].x, A[18].x, b_l.x + b.x
    else:
        return False, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None


# min hNN-gNN
def relax_verify_initCond(wb_1_list, wb_2_list, wb_length, model):
    cex_init = []
    res_init = []

    init_len = 0.3
    cnt = 0
    cex_idx_list = []

    while True:
        relax_flag = False

        new_cex_list = []
        if len(cex_idx_list) == 0:
            cex_idx_list.append([init_min, init_max])

        # print("init_cex_idx_list_num:", len(cex_idx_list))

        for cex_idx in cex_idx_list:
            block_min, block_max = cex_idx
            x0_l = block_min[0]
            x1_l = block_min[1]
            x2_l = block_min[2]
            x3_l = block_min[3]
            x4_l = block_min[4]
            x5_l = block_min[5]
            x6_l = block_min[6]
            x7_l = block_min[7]
            x8_l = block_min[8]
            x9_l = block_min[9]
            x10_l = block_min[10]
            x11_l = block_min[11]
            x12_l = block_min[12]
            x13_l = block_min[13]
            x14_l = block_min[14]
            x15_l = block_min[15]
            x16_l = block_min[16]
            x17_l = block_min[17]
            x18_l = block_min[18]
            x0_r = block_max[0]
            x1_r = block_max[1]
            x2_r = block_max[2]
            x3_r = block_max[3]
            x4_r = block_max[4]
            x5_r = block_max[5]
            x6_r = block_max[6]
            x7_r = block_max[7]
            x8_r = block_max[8]
            x9_r = block_max[9]
            x10_r = block_max[10]
            x11_r = block_max[11]
            x12_r = block_max[12]
            x13_r = block_max[13]
            x14_r = block_max[14]
            x15_r = block_max[15]
            x16_r = block_max[16]
            x17_r = block_max[17]
            x18_r = block_max[18]

            m = gp.Model()
            m.setParam('Outputflag', 0)
            x0 = m.addMVar((DIM,), vtype=GRB.CONTINUOUS, lb=-np.inf, ub=np.inf)

            m.addConstr(x0[0] >= x0_l)
            m.addConstr(x0[0] <= x0_r)
            m.addConstr(x0[1] >= x1_l)
            m.addConstr(x0[1] <= x1_r)
            m.addConstr(x0[2] >= x2_l)
            m.addConstr(x0[2] <= x2_r)
            m.addConstr(x0[3] >= x3_l)
            m.addConstr(x0[3] <= x3_r)
            m.addConstr(x0[4] >= x4_l)
            m.addConstr(x0[4] <= x4_r)
            m.addConstr(x0[5] >= x5_l)
            m.addConstr(x0[5] <= x5_r)
            m.addConstr(x0[6] >= x6_l)
            m.addConstr(x0[6] <= x6_r)
            m.addConstr(x0[7] >= x7_l)
            m.addConstr(x0[7] <= x7_r)
            m.addConstr(x0[8] >= x8_l)
            m.addConstr(x0[8] <= x8_r)
            m.addConstr(x0[9] >= x9_l)
            m.addConstr(x0[9] <= x9_r)
            m.addConstr(x0[10] >= x10_l)
            m.addConstr(x0[10] <= x10_r)
            m.addConstr(x0[11] >= x11_l)
            m.addConstr(x0[11] <= x11_r)
            m.addConstr(x0[12] >= x12_l)
            m.addConstr(x0[12] <= x12_r)
            m.addConstr(x0[13] >= x13_l)
            m.addConstr(x0[13] <= x13_r)
            m.addConstr(x0[14] >= x14_l)
            m.addConstr(x0[14] <= x14_r)
            m.addConstr(x0[15] >= x15_l)
            m.addConstr(x0[15] <= x15_r)
            m.addConstr(x0[16] >= x16_l)
            m.addConstr(x0[16] <= x16_r)
            m.addConstr(x0[17] >= x17_l)
            m.addConstr(x0[17] <= x17_r)
            m.addConstr(x0[18] >= x18_l)
            m.addConstr(x0[18] <= x18_r)

            hNN = barrier_f_without_der(m, x0, wb_2_list, wb_length, True)

            # hat_gNN = Ax + b = hNN upper bound
            feasible, A0, A1, A2, A3, A4, A5, A6, A7, A8, A9, A10, A11, A12, A13, A14, A15, A16, A17, A18, b = barrier_f_linear(wb_1_list, wb_length, [x0_l, x1_l, x2_l, x3_l, x4_l, x5_l, x6_l, x7_l, x8_l, x9_l, x10_l, x11_l, x12_l, x13_l, x14_l, x15_l, x16_l, x17_l, x18_l], [x0_r, x1_r, x2_r, x3_r, x4_r, x5_r, x6_r, x7_r, x8_r, x9_r, x10_r, x11_r, x12_r, x13_r, x14_r, x15_r, x16_r, x17_r, x18_r])
            if not feasible:
                # print("unfeasible init barrier_f_linear")
                continue

            m.setObjective(hNN - (A0 * x0[0] + A1 * x0[1] + A2 * x0[2] + A3 * x0[3] + A4 * x0[4] + A5 * x0[5] + A6 * x0[6] + A7 * x0[7] + A8 * x0[8] + A9 * x0[9] + A10 * x0[10]
                                  + A11 * x0[11] + A12 * x0[12] + A13 * x0[13] + A14 * x0[14] + A15 * x0[15] + A16 * x0[16] + A17 * x0[17] + A18 * x0[18] + b), GRB.MINIMIZE)
            m.optimize()

            if m.status == 2:
                if m.objVal < 0:
                    if model(torch.tensor(np.array(x0.x), dtype=torch.float64)) > 0:
                        cex_init.append(np.array(x0.x))
                        res_init.append(m.objVal)
                        break
                    else:
                        relax_flag = True
                        new_cex_list += block_space(block_min, block_max)

        if relax_flag and len(cex_init) == 0 and cnt < cnt_num:
            init_len /= 2
            cnt += 1
            print("init_len:", init_len)
            cex_idx_list = new_cex_list

        else:
            break

    return cex_init, res_init

# min gNN-hNN
def relax_verify_unsafeCond(wb_1_list, wb_2_list, wb_length, model):
    cex_unsafe = []
    res_unsafe = []

    unsafe_len = 0.05
    cnt = 0
    cex_idx_list = []
    while True:
        relax_flag = False

        new_cex_list = []
        if len(cex_idx_list) == 0:
            cex_idx_list.append([unsafe_min, unsafe_max])

        # print("unsafe_cex_idx_list_num:", len(cex_idx_list))

        for cex_idx in cex_idx_list:
            block_min, block_max = cex_idx
            x0_l = block_min[0]
            x1_l = block_min[1]
            x2_l = block_min[2]
            x3_l = block_min[3]
            x4_l = block_min[4]
            x5_l = block_min[5]
            x6_l = block_min[6]
            x7_l = block_min[7]
            x8_l = block_min[8]
            x9_l = block_min[9]
            x10_l = block_min[10]
            x11_l = block_min[11]
            x12_l = block_min[12]
            x13_l = block_min[13]
            x14_l = block_min[14]
            x15_l = block_min[15]
            x16_l = block_min[16]
            x17_l = block_min[17]
            x18_l = block_min[18]
            x0_r = block_max[0]
            x1_r = block_max[1]
            x2_r = block_max[2]
            x3_r = block_max[3]
            x4_r = block_max[4]
            x5_r = block_max[5]
            x6_r = block_max[6]
            x7_r = block_max[7]
            x8_r = block_max[8]
            x9_r = block_max[9]
            x10_r = block_max[10]
            x11_r = block_max[11]
            x12_r = block_max[12]
            x13_r = block_max[13]
            x14_r = block_max[14]
            x15_r = block_max[15]
            x16_r = block_max[16]
            x17_r = block_max[17]
            x18_r = block_max[18]

            m = gp.Model()
            m.setParam('Outputflag', 0)
            x0 = m.addMVar((DIM,), vtype=GRB.CONTINUOUS, lb=-np.inf, ub=np.inf)

            m.addConstr(x0[0] >= x0_l)
            m.addConstr(x0[0] <= x0_r)
            m.addConstr(x0[1] >= x1_l)
            m.addConstr(x0[1] <= x1_r)
            m.addConstr(x0[2] >= x2_l)
            m.addConstr(x0[2] <= x2_r)
            m.addConstr(x0[3] >= x3_l)
            m.addConstr(x0[3] <= x3_r)
            m.addConstr(x0[4] >= x4_l)
            m.addConstr(x0[4] <= x4_r)
            m.addConstr(x0[5] >= x5_l)
            m.addConstr(x0[5] <= x5_r)
            m.addConstr(x0[6] >= x6_l)
            m.addConstr(x0[6] <= x6_r)
            m.addConstr(x0[7] >= x7_l)
            m.addConstr(x0[7] <= x7_r)
            m.addConstr(x0[8] >= x8_l)
            m.addConstr(x0[8] <= x8_r)
            m.addConstr(x0[9] >= x9_l)
            m.addConstr(x0[9] <= x9_r)
            m.addConstr(x0[10] >= x10_l)
            m.addConstr(x0[10] <= x10_r)
            m.addConstr(x0[11] >= x11_l)
            m.addConstr(x0[11] <= x11_r)
            m.addConstr(x0[12] >= x12_l)
            m.addConstr(x0[12] <= x12_r)
            m.addConstr(x0[13] >= x13_l)
            m.addConstr(x0[13] <= x13_r)
            m.addConstr(x0[14] >= x14_l)
            m.addConstr(x0[14] <= x14_r)
            m.addConstr(x0[15] >= x15_l)
            m.addConstr(x0[15] <= x15_r)
            m.addConstr(x0[16] >= x16_l)
            m.addConstr(x0[16] <= x16_r)
            m.addConstr(x0[17] >= x17_l)
            m.addConstr(x0[17] <= x17_r)
            m.addConstr(x0[18] >= x18_l)
            m.addConstr(x0[18] <= x18_r)

            gNN = barrier_f_without_der(m, x0, wb_1_list, wb_length, True)

            # hat_hNN = Ax + b = hNN upper bound
            feasible, A0, A1, A2, A3, A4, A5, A6, A7, A8, A9, A10, A11, A12, A13, A14, A15, A16, A17, A18, b = barrier_f_linear(wb_2_list, wb_length, [x0_l, x1_l, x2_l, x3_l, x4_l, x5_l, x6_l, x7_l, x8_l, x9_l, x10_l, x11_l, x12_l, x13_l, x14_l, x15_l, x16_l, x17_l, x18_l], [x0_r, x1_r, x2_r, x3_r, x4_r, x5_r, x6_r, x7_r, x8_r, x9_r, x10_r, x11_r, x12_r, x13_r, x14_r, x15_r, x16_r, x17_r, x18_r])
            if not feasible:
                # print("unfeasible unsafe barrier_f_linear")
                continue

            m.setObjective(gNN - (A0 * x0[0] + A1 * x0[1] + A2 * x0[2] + A3 * x0[3] + A4 * x0[4] + A5 * x0[5] + A6 * x0[6] + A7 * x0[7] + A8 * x0[8] + A9 * x0[9] + A10 * x0[10]
                                  + A11 * x0[11] + A12 * x0[12] + A13 * x0[13] + A14 * x0[14] + A15 * x0[15] + A16 * x0[16] + A17 * x0[17] + A18 * x0[18] + b), GRB.MINIMIZE)
            m.optimize()

            # print("minimize unsafe value:", m.objVal)
            if m.status == 2:
                if m.objVal < 0:
                    if model(torch.tensor(np.array(x0.x), dtype=torch.float64)) < 0:
                        cex_unsafe.append(np.array(x0.x))
                        res_unsafe.append(m.objVal)
                        break
                    else:
                        relax_flag = True
                        new_cex_list += block_space(block_min, block_max)

        if relax_flag and len(cex_unsafe) == 0 and cnt < cnt_num:
            unsafe_len /= 2
            cnt += 1
            print("unsafe_len:", unsafe_len)
            cex_idx_list = new_cex_list

        else:
            break
    return cex_unsafe, res_unsafe


# max ( (gNN(f)-hNN(f)) - (1-lambda) * (gNN-hNN) )
def relax_verify_thirdCond(wb_1_list, wb_2_list, wb_length, x_len, model):
    cex_domain = []
    res_domain = []

    domain_len = 0.6
    cnt = 0

    x0_min = domain_min[0]
    x0_max = domain_max[0]
    x1_min = domain_min[1]
    x1_max = domain_max[1]
    x2_min = domain_min[2]
    x2_max = domain_max[2]
    x3_min = domain_min[3]
    x3_max = domain_max[3]
    x4_min = domain_min[4]
    x4_max = domain_max[4]
    x5_min = domain_min[5]
    x5_max = domain_max[5]
    x6_min = domain_min[6]
    x6_max = domain_max[6]
    x7_min = domain_min[7]
    x7_max = domain_max[7]
    x8_min = domain_min[8]
    x8_max = domain_max[8]
    x9_min = domain_min[9]
    x9_max = domain_max[9]
    x10_min = domain_min[10]
    x10_max = domain_max[10]
    x11_min = domain_min[11]
    x11_max = domain_max[11]
    x12_min = domain_min[12]
    x12_max = domain_max[12]
    x13_min = domain_min[13]
    x13_max = domain_max[13]
    x14_min = domain_min[14]
    x14_max = domain_max[14]
    x15_min = domain_min[15]
    x15_max = domain_max[15]
    x16_min = domain_min[16]
    x16_max = domain_max[16]
    x17_min = domain_min[17]
    x17_max = domain_max[17]
    x18_min = domain_min[18]
    x18_max = domain_max[18]

    cex_idx_list = []
    while True:
        relax_flag = False

        new_cex_list = []
        if len(cex_idx_list) == 0:
            cex_idx_list.append([domain_min, domain_max])

        print("domain_cex_idx_list_num:", len(cex_idx_list))

        for cex_idx in cex_idx_list:
            block_min, block_max = cex_idx
            x0_l = block_min[0]
            x1_l = block_min[1]
            x2_l = block_min[2]
            x3_l = block_min[3]
            x4_l = block_min[4]
            x5_l = block_min[5]
            x6_l = block_min[6]
            x7_l = block_min[7]
            x8_l = block_min[8]
            x9_l = block_min[9]
            x10_l = block_min[10]
            x11_l = block_min[11]
            x12_l = block_min[12]
            x13_l = block_min[13]
            x14_l = block_min[14]
            x15_l = block_min[15]
            x16_l = block_min[16]
            x17_l = block_min[17]
            x18_l = block_min[18]
            x0_r = block_max[0]
            x1_r = block_max[1]
            x2_r = block_max[2]
            x3_r = block_max[3]
            x4_r = block_max[4]
            x5_r = block_max[5]
            x6_r = block_max[6]
            x7_r = block_max[7]
            x8_r = block_max[8]
            x9_r = block_max[9]
            x10_r = block_max[10]
            x11_r = block_max[11]
            x12_r = block_max[12]
            x13_r = block_max[13]
            x14_r = block_max[14]
            x15_r = block_max[15]
            x16_r = block_max[16]
            x17_r = block_max[17]
            x18_r = block_max[18]

            m = gp.Model()
            m.setParam('Outputflag', 0)
            m.setParam('NonConvex', 2)
            # m.setParam('InfUnbdInfo', 1)

            x0 = m.addMVar((DIM,), vtype=GRB.CONTINUOUS, name='x0', lb=-np.inf, ub=np.inf)
            x0[0].lb = x0_l
            x0[0].ub = x0_r
            x0[1].lb = x1_l
            x0[1].ub = x1_r
            x0[2].lb = x2_l
            x0[2].ub = x2_r
            x0[3].lb = x3_l
            x0[3].ub = x3_r
            x0[4].lb = x4_l
            x0[4].ub = x4_r
            x0[5].lb = x5_l
            x0[5].ub = x5_r
            x0[6].lb = x6_l
            x0[6].ub = x6_r
            x0[7].lb = x7_l
            x0[7].ub = x7_r
            x0[8].lb = x8_l
            x0[8].ub = x8_r
            x0[9].lb = x9_l
            x0[9].ub = x9_r
            x0[10].lb = x10_l
            x0[10].ub = x10_r
            x0[11].lb = x11_l
            x0[11].ub = x11_r
            x0[12].lb = x12_l
            x0[12].ub = x12_r
            x0[13].lb = x13_l
            x0[13].ub = x13_r
            x0[14].lb = x14_l
            x0[14].ub = x14_r
            x0[15].lb = x15_l
            x0[15].ub = x15_r
            x0[16].lb = x16_l
            x0[16].ub = x16_r
            x0[17].lb = x17_l
            x0[17].ub = x17_r
            x0[18].lb = x18_l
            x0[18].ub = x18_r

            x1 = m.addMVar((DIM,), vtype=GRB.CONTINUOUS, name='x1', lb=-np.inf, ub=np.inf)  # x1 = F(kx+b; x0)
            x1[0].lb = x0_min
            x1[0].ub = x0_max
            x1[1].lb = x1_min
            x1[1].ub = x1_max
            x1[2].lb = x2_min
            x1[2].ub = x2_max
            x1[3].lb = x3_min
            x1[3].ub = x3_max
            x1[4].lb = x4_min
            x1[4].ub = x4_max
            x1[5].lb = x5_min
            x1[5].ub = x5_max
            x1[6].lb = x6_min
            x1[6].ub = x6_max
            x1[7].lb = x7_min
            x1[7].ub = x7_max
            x1[8].lb = x8_min
            x1[8].ub = x8_max
            x1[9].lb = x9_min
            x1[9].ub = x9_max
            x1[10].lb = x10_min
            x1[10].ub = x10_max
            x1[11].lb = x11_min
            x1[11].ub = x11_max
            x1[12].lb = x12_min
            x1[12].ub = x12_max
            x1[13].lb = x13_min
            x1[13].ub = x13_max
            x1[14].lb = x14_min
            x1[14].ub = x14_max
            x1[15].lb = x15_min
            x1[15].ub = x15_max
            x1[16].lb = x16_min
            x1[16].ub = x16_max
            x1[17].lb = x17_min
            x1[17].ub = x17_max
            x1[18].lb = x18_min
            x1[18].ub = x18_max

            xx1 = m.addVar()
            xx3 = m.addVar()
            xx5 = m.addVar()
            xx7 = m.addVar()
            xx9 = m.addVar()
            xx11 = m.addVar()
            xx13 = m.addVar()
            xx15 = m.addVar()
            xx17 = m.addVar()
            m.addConstr(xx1 == x0[1] * x0[1])
            m.addConstr(xx3 == x0[3] * x0[3])
            m.addConstr(xx5 == x0[5] * x0[5])
            m.addConstr(xx7 == x0[7] * x0[7])
            m.addConstr(xx9 == x0[9] * x0[9])
            m.addConstr(xx11 == x0[11] * x0[11])
            m.addConstr(xx13 == x0[13] * x0[13])
            m.addConstr(xx15 == x0[15] * x0[15])
            m.addConstr(xx17 == x0[17] * x0[17])
            m.addConstr(x1[0] == x0[0] + time_step * ((x0[1] + 2 * x0[2] + 2 * x0[3] + 2 * x0[4] + x0[5]) / 100.0 + 1))
            m.addConstr(x1[1] == x0[1] + time_step * (x0[2]))
            m.addConstr(x1[2] == x0[2] + time_step * (-10 * (x0[1] - (xx1 * x0[1]) / 6.0) - x0[1]))
            m.addConstr(x1[3] == x0[3] + time_step * (x0[4]))
            m.addConstr(x1[4] == x0[4] + time_step * (-10 * (x0[3] - (xx3 * x0[3]) / 6.0) - x0[1]))
            m.addConstr(x1[5] == x0[5] + time_step * (x0[6]))
            m.addConstr(x1[6] == x0[6] + time_step * (-10 * (x0[5] - (xx5 * x0[5]) / 6.0) - x0[1]))
            m.addConstr(x1[7] == x0[7] + time_step * (x0[8]))
            m.addConstr(x1[8] == x0[8] + time_step * (-10 * (x0[7] - (xx7 * x0[7]) / 6.0) - x0[1]))
            m.addConstr(x1[9] == x0[9] + time_step * (x0[10]))
            m.addConstr(x1[10] == x0[10] + time_step * (-10 * (x0[9] - (xx9 * x0[9]) / 6.0) - x0[1]))
            m.addConstr(x1[11] == x0[11] + time_step * (x0[12]))
            m.addConstr(x1[12] == x0[12] + time_step * (-10 * (x0[11] - (xx11 * x0[11]) / 6.0) - x0[1]))
            m.addConstr(x1[13] == x0[13] + time_step * (x0[14]))
            m.addConstr(x1[14] == x0[14] + time_step * (-10 * (x0[13] - (xx13 * x0[13]) / 6.0) - x0[1]))
            m.addConstr(x1[15] == x0[15] + time_step * (x0[16]))
            m.addConstr(x1[16] == x0[16] + time_step * (-10 * (x0[15] - (xx15 * x0[15]) / 6.0) - x0[1]))
            m.addConstr(x1[17] == x0[17] + time_step * (x0[18]))
            m.addConstr(x1[18] == x0[18] + time_step * (-10 * (x0[17] - (xx17 * x0[17]) / 6.0) - x0[1]))

            # get_barrier_start = time.time()
            gNN = barrier_f_without_der(m, x0, wb_1_list, wb_length, True)
            hNN_f = barrier_f_without_der(m, x1, wb_2_list, wb_length, True)
            # get_barrier_end = time.time()
            # print("get_barrier：", get_barrier_end - get_barrier_start)

            # get_linear_start = time.time()
            feasible, gNN_f_A0, gNN_f_A1, gNN_f_A2, gNN_f_A3, gNN_f_A4, gNN_f_A5, gNN_f_A6, gNN_f_A7, gNN_f_A8, gNN_f_A9, gNN_f_A10, gNN_f_A11, gNN_f_A12, gNN_f_A13, gNN_f_A14, gNN_f_A15, gNN_f_A16, gNN_f_A17, gNN_f_A18, gNN_b = barrier_Bf_nonlinear(wb_1_list, wb_length, [x0_l, x1_l, x2_l, x3_l, x4_l, x5_l, x6_l, x7_l, x8_l, x9_l, x10_l, x11_l, x12_l, x13_l, x14_l, x15_l, x16_l, x17_l, x18_l], [x0_r, x1_r, x2_r, x3_r, x4_r, x5_r, x6_r, x7_r, x8_r, x9_r, x10_r, x11_r, x12_r, x13_r, x14_r, x15_r, x16_r, x17_r, x18_r])
            if not feasible:
                continue
            feasible, hNN_A0, hNN_A1, hNN_A2, hNN_A3, hNN_A4, hNN_A5, hNN_A6, hNN_A7, hNN_A8, hNN_A9, hNN_A10, hNN_A11, hNN_A12, hNN_A13, hNN_A14, hNN_A15, hNN_A16, hNN_A17, hNN_A18, hNN_b = barrier_f_linear(wb_2_list, wb_length, [x0_l, x1_l, x2_l, x3_l, x4_l, x5_l, x6_l, x7_l, x8_l, x9_l, x10_l, x11_l, x12_l, x13_l, x14_l, x15_l, x16_l, x17_l, x18_l], [x0_r, x1_r, x2_r, x3_r, x4_r, x5_r, x6_r, x7_r, x8_r, x9_r, x10_r, x11_r, x12_r, x13_r, x14_r, x15_r, x16_r, x17_r, x18_r])
            if not feasible:
                continue
            # get_linear_end = time.time()
            # print("get_linear：", get_linear_end - get_linear_start)

            m.setObjective( (gNN_f_A0*x1[0] + gNN_f_A1*x1[1] + gNN_f_A2*x1[2] + gNN_f_A3*x1[3] + gNN_f_A4*x1[4] + gNN_f_A5*x1[5] + gNN_f_A6*x1[6] + gNN_f_A7*x1[7] + gNN_f_A8*x1[8] + gNN_f_A9*x1[9]
                             + gNN_f_A10*x1[10] + gNN_f_A11*x1[11] + gNN_f_A12*x1[12] + gNN_f_A13*x1[13] + gNN_f_A14*x1[14] + gNN_f_A15*x1[15] + gNN_f_A16*x1[16] + gNN_f_A17*x1[17] + gNN_f_A18*x1[18] + gNN_b)
                            -hNN_f
                            - (1-alpha) * gNN
                            + (1-alpha) * (hNN_A0*x0[0] + hNN_A1*x0[1] + hNN_A2*x0[2] + hNN_A3*x0[3] + hNN_A4*x0[4] + hNN_A5*x0[5] + hNN_A6*x0[6] + hNN_A7*x0[7] + hNN_A8*x0[8] + hNN_A9*x0[9]
                                           + hNN_A10*x0[10] + hNN_A11*x0[11] + hNN_A12*x0[12] + hNN_A13*x0[13] + hNN_A14*x0[14] + hNN_A15*x0[15] + hNN_A16*x0[16] + hNN_A17*x0[17] + hNN_A18*x0[18] + hNN_b),
                            GRB.MAXIMIZE)

            m.optimize()
            # print(m.status)
            if m.status == 2:
                # print(m.objVal)
                if m.objVal > 0:
                    # print(f"m.objVal:{m.objVal}, model(x0.x):{model(vector_field(torch.tensor([np.array(x0.x)]))) - (1-alpha) * model(torch.tensor(np.array(x0.x)))}")
                    if model(vector_field(torch.tensor(np.array(x0.x)[None, :]))) - (1 - alpha) * model(
                            torch.tensor(np.array(x0.x))) > 0:
                        cex_domain.append(np.array(x0.x))
                        res_domain.append(m.objVal)
                        # print(np.array(x0.x))
                        # print(cex_domain)
                        # print(m.objVal)
                        # print(res_domain)
                        break
                    else:
                        relax_flag = True
                        if cnt < 1:
                            new_cex_list += block_space(block_min, block_max)

        if relax_flag and len(cex_domain) == 0 and cnt < cnt_num:
            domain_len /= 2
            cnt += 1
            print("domain_len:", domain_len)
            cex_idx_list = new_cex_list
        else:
            break

    return cex_domain, res_domain


def verif(model, x_len):
    wb_1_list = []
    wb_2_list = []
    params = model.state_dict()
    i = 0
    for k in params:
        # print(k)
        # 'hidden_layers1.0.weight' 'hidden_layers1.0.bias' 'hidden_layers2.0.weight' 'hidden_layers2.0.bias'
        if 'layers' in k:
            continue
        if i % 4 < 2:
            wb_1_list.append(params[k].detach().numpy())
        else:
            wb_2_list.append(params[k].detach().numpy())
        i += 1

    wb_length = int(len(wb_1_list) / 2)
    # time_init_start = time.time()
    cex_init, res_init = relax_verify_initCond(wb_1_list, wb_2_list, wb_length, model)
    # time_init_end = time.time()
    # print("init cond time:", time_init_end - time_init_start)

    # time_unsafe_start = time.time()
    cex_unsafe, res_unsafe = relax_verify_unsafeCond(wb_1_list, wb_2_list, wb_length, model)
    # time_unsafe_end = time.time()
    # print("unsafe cond time:", time_unsafe_end - time_unsafe_start)

    # time_domain_start = time.time()
    cex_domain, res_domain = relax_verify_thirdCond(wb_1_list, wb_2_list, wb_length, x_len, model)
    # time_domain_end = time.time()
    # print("domain cond time:", time_domain_end - time_domain_start)

    print("cex_init: ", cex_init)
    print("res_init: ", res_init)
    print("res_unsafe: ", res_unsafe)
    print("cex_unsafe: ", cex_unsafe)
    print("cex_domain: ", cex_domain)
    print("res_domain: ", res_domain)
    return cex_init, cex_unsafe, cex_domain

